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A  computer  program  for  the  calculation  of  linearised  steady  and 
oscillatory  supersonic  flow  over  thin  wings  has  been  implemented  on 
the  ELXSI.  The  program  uses  an  explicit  finite  difference  scheme  on 
a  characteristic  grid  to  solve  for  the  reduced  potential  and  pressure 
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NOTATION 


Symbol 

Definition 

AME 

final  Mach  number 

AMS 

initial  Mach  number 

AOME 

final  reduced  frequency 

AOMS 

initial  reduced  frequency 

BZ 

semi  span  length 

C(i,J  ,k) 

f(i.J-1.  k-1)  +  k+1)  ^  f(i,J  +  1,  k-1) 

+  f(l.J  +  l .  k+1 ) 

S 

steady  pressure  coefficient 

DAM 

delta  change  in  Mach  number 

DAOM 

delta  change  in  reduced  frequency 

DX 

characteristic  grid  spacing 

f 

reduced  steady/oscillatory  velocity  potential 

f(i.j,l<) 

(1-p)  f(l,J,k)  +  1  [f(i+l.  J,  k)  ♦  f{i-1,  J 

FXL  (Y) 

position  of  the  leading  edge 

FXT  (Y) 

position  of  the  trailing  edge 

h 

vertical  position  of  wing  surface 

f’o 

oscillatory  part  of  wing  surface 

steady  upper /lower  position  of  wing 

HI  (X,Y) 

hg/hg*j5  for  oscillatory/steady  flow 

H2  (X,T) 

hox/Hs”j{  for  oscillatory/steady  flow 

lOUT 

output  parameter 

ISW 

option  parameter 

k 

HM/o^ 

K 

k/zero  for  oscillatory/steady  flow 

1 


typical  chord  length 


NOTATION  (CONT.) 


Symbol 


Definition 


M 

M(i.j  ,k) 


P 

P 

a 

Q 


U 

CD 


x.y.z.t 


a 

6 

A 


n,;,T 

X 

p« 


« 

I- 

u 


free  stream  Mach  number 

f(i.j-l.k)  +  f(i,j*1,k)  +  f(i,j,k-l) 

+  f(l,J.k+1) 

constant/local  pressure 

free  stream  pressure 

oscillatory  pressure  coefficient 

free  stream  velocity 

Euler ian  variables 

/  -1 

thickness  ratio 
rectangular  grid  step  length 
Prandtl-Glauert  variables 
constant 

free  stream  density 

perturbation  potential 

oscillatory  perturbation  potential 

steady  perturbation  potential 

sealed  oscillatory  perturbation  potential 

velocity  potential 


reduced  frequency 


-IP-- 


▼ 


(1) 


1 .  INTRODUCTION 


A  computer  program  for  the  calculation  of  steady  and  oscillatory 
supersonic  flow  over  thin  wings  has  been  implemented  in  Fortran  on  the 
ELXSI.  As  the  basic  differential  equation  is  'time-like'  and  the  relevant 
part  of  the  flow-field  finite,  the  program  uses  an  explicit  finite- 
difference  scheme  to  solve  for  the  reduced  potential  on  the  wing. 


In  section  2  the  governing  differential  equations  and  boundary 
conditions  for  steady  and  oscillatory  supersonic  flow  over  a  thin  wing  are 
developed,  and  in  section  3  the  explicit  finite  difference  formula 
applicable  to  linearised  supersonic  flow  is  presented.  This  formula  has 
the  advantage  of  using  a  characteristic  grid  rather  than  a  rectangular 
grid  thus  halving  the  number  of  grid  points.  The  method  is  restricted  to 
boundary  conditions  defined  on  a  plane.  Section  it  discusses  how  to  use 
the  program,  that  is;  setting  up  the  data,  compiling,  binding  and 
running,  A  listing  of  the  program  is  given  in  Appendix  A  and  finally  a 
listing  of  a  subroutine  to  read  the  output  from  the  supersonic  program  is 
given  in  Appendix  B,  this  subroutine  is  discussed  in  Section  it. 


2.  GOVERNING  EQUATION  AND  BOUNDART  CONDITIONS 


In  the  following,  the  isentropic  invisoid  flow  of  a  perfect  gas, 
initially  irrotational ,  is  considered.  Thus  we  may  assume  the  ejcistence 
of  a  velocity  potential  i(i,  such  that  the  fluid  velocity  v  -  yi(i.  If 
U  Is  the  free  stream  velocity  (which  is  assumed  to  be  directed  solely 
in  the  X-  direction),  then  a  perturbation  potential  it>  can  be  defined 
such  that 


(x,y,z,t)  =  U~  [x  +  $  (x,y,z,t)],  (1) 


where  x,y,z  represent  a  rectangular  carterlan  co-ordinate  system  and  t  is 
the  time  variable.  Note  that  subsequent  equations  are  expressed  in  non- 
dimensional  co-ordinates  based  on  a  length  scale  I,  a  typical  value  of 
the  airfoil  chord  length,  a  velocity  scale  U  and  a  density  scale 
a  typical  value  of  the  density  at  infinity.  "-The  time  variable  is  then 
scaled  with  1/U^  and  the  pressure  with 


For  thin  nearly  planar  wings  of  moderate  aspect  ratio  the  governing 
differential  equation  for  the  perturbation  velocity  potential  can  be 
written  as 


♦xx  "  Vy  "  *22  *  '^xt  *  ■  °' 


(2) 


I 


(2) 


T 


where  M  is  the  free  stream  Mach  number.  It  should  be  noted  here  that  the 
procedure  used  to  derive  (2)  is  based  upon  the  small  parameter  6 
representing  the  ratio  of  airfoil  thickness  to  chord  length.  It  is 
assumed  that  as  S  -*  o,  (ti/6  remains  fixed.  Then  ignoring  terms  of 

second  order  in  6,  we  find  that  (2)  is  the  appropriate  governing 

differential  equation  for 


On  a  surface  in  an  inviscid  fluid  the  usual  kinematic  condition 
applies.  If  Z  =  h(x,y,t)  specifies  the  wing  (where  h  is  under  6  by 
definition),  then,  to  first  order  in  6,  the  boundary  condition  on  the 
wing  is 


z=0 


3x  *  at  * 


(3) 


In  many  practical  applications  the  unsteady  motion  of  a  wing  may  be 
assumed  to  consist  of  small  inf initerimal  perturbations  around  the  steady 
state  configuration.  Thus  the  motion  consists  of  a  steady  component  plus 
a  small  harmonically  oscillating  unsteady  component.  Let 


«  (x,y,z,t)  -  (Jig  (x,y,z)  +  <|.^  (x,y,z) 

('la) 

h  (x,y,t)  -  hg"  (x,y)  +  h^  (x,y) 

(Hb) 

where,  the  subscript  S/0  denotes  the  steady/osoillatory  part,  h  ” 
represents  the  upper /lower  surface  of  the  wing  and  due  to  the  sealings 
introduced  above,  u  represents  the  reduced  frequency  or  Stroubal 

number.  Note  that  the  actual  frequency  is  u  /  ,.  Substituting  (ka 

and  b)  into  (2)  and  (3),  and  separating  steady  and  ^oscillatory  parts  we 
find  that 

2 

^  *sxx  '  '*'syy  *  '*’szz’ 

(5a) 

♦sz  ^  Z.0  ■  h 

(5b) 

•  • 

and 

’  ♦oxx  *  ^i“  ^ox  ■  *0 

'^oyy  *  ^ozz’ 

(6a) 

3,. 

3ho 

*oz  ^  z-o  ■  aT  *  ^ 

(6b) 

1 


(3) 


At  this  stage  it  is  convenient  to  introduce  the  generalised  Prandtl- 
Glauert  variables 

2 

n  =  ay,  ;  =  az  and  x  =  a  t/m  (7a) 

where  a^  =  -=1  >0.  (7b) 


Also,  in  order  to  write  (6a)  in  canonical  form  let 


=  ♦(x.n.c) 


-ikmx 


(8a) 


where 


k 


u  M/o 


2 


(8b) 


Substituting  (7a  and  b)  into  (5a  and  b),  and,  (7a  and  b)  and  (8a  and  b) 
into  (6a  and  b)  the  steady  and  oscillatory  problems  become 


A 

^SXX 


’^snn 


(9a) 


'•'sc/;  -  0± 


1 


a 


on  the  wing 


(9b) 


and 


hh 


(10a) 


ikmx 

V?-0  ■  f  ^^x  *  ^“ho^ 

Finally  note  that  if  k  is  zero  then  the  partial  differential  equation 
(10a)  reduces  to  the  aai,.j  form  as  (9a).  Also  the  boundary  conditions  (9b) 
and  (10b)  have  similar  forms.  Due  to  these  similarities  both  the  steady 
and  oscillatory  problems  can  be  solved  by  the  same  finite  difference 
technique.  In  the  next  section  we  discuss  the  finite  difference  scheme 
that  will  be  used  to  solve  both  (9a, b)  and  (10a, b). 


iL 


3.  FINITE  DIFFERENCE  SCHEME 


The  basic  equation  of  linearised  supersonic  flow  is 


f  +  K  f 

XX 


f  f 

yy  zz 


where  x  is  the  flow  direction,  f  the  reduced  perturbation  velocity 
potential  and  K  is  a  frequency  parameter  for  unsteady  flow,  or  zero  for 
steady  flow. 


Let  the  potential  be  defined  at  grid  points  forming  a  cubic  lattice 
at  which  it  may  be  assumed  that  the  co-ordinates  take  integer  values 
(i,J,k).  It  has  been  shown  by  Sullivan  [1]  that  a  second  order, 

consistent,  finite  difference  scheme,  with  rotational  symmetry  about  the 
x-axis  and  which  reduces  to  the  method  of  characteristics  for  two- 
dimensional  flow  in  the  x-y  or  x-z  planes,  is 


f(i+1.J,k)  +  f(l-1,J,k)  ♦  f  (i,J.k) 

-  -1  f(i,J,k)  +  4  M(i,j,k)  +(4-4)  C(i,j,k)  (12a) 


f(i.J.k)  -  (1-p)  f(l,J.k)  +  4  P  Cf(i  +  1,J,k)  +  f(i-1,J,k)],  (12b) 


M(i,J,k)  -  f(i,j-l,k)  +  f(i,j  +  1,k)  +  f(l,j,k-l) 

+  f(l,j,K+1). 

C(l,J,k)  =.  f(i,J-l,k-1)  +  f(i,J+1,k-l) 

+  f(i,j-1,k+l)  ♦  f(i,J+1,k+l), 


where  P  and  1  are  constants  and  A  is  the  grid  step-length.  It  has 
also  been  shown  by  Sullivan  [1]  that  the  two-parameter  family  of  finite 
difference  schemes  (12a, b,c  and  d)  are  strictly  stable  in  the  range 


12  10 


Pi  2 


giving: 


Following  Sullivan  [1J  choose  the  simplest  scheme  1-0,  p  -  1 


1 


tlii) 


▼ 


T 


f(Ul  .j  .k)  +  f(i-l  .j  .k) 


(2»K^A‘^) 


C(i  ,J  ,k) . 


Note  that  (1^)  iias  the  advantage  that  the  terras  f{i,J,k)  and  M(i,j,k)  do 
not  appear,  giving  a  reduction  in  the  number  of  grid  points  required. 
That  is,  the  system  (l*t)  is  solved  on  a  characteristic  grid  rather  than  a 
rectangular  grid. 


Note  that  (1^)  is  neutrally  stable.  This  neutral  stability  will  raanifest 
itself  as  numerical  instability  as  the  number  of  grid  points  downstream 
becomes  large. 


Normal  derivative  boundary  conditions  can  be  implemented  in  a 
variety  of  ways.  The  technique  used  with  (14)  is  to  generate  dummy 
potentials  using  the  known  normal  derivates,  with  grid  points  on  the 
boundary  treated  as  interior  points.  That  is 


fz  /  =  [f(i,j.2)  -  f(i.j,0)]  /  2A, 


(15) 


then  on  the  boundary,  level  k  =  1,  (14)  becomes 


f(iM,j,l) 

2 

“  - FF 

(2+K  A  ) 


f(l-1  ,j  ,1  ) 

f(i.J-l,2)  ♦  f(i,j*1,2)  -  A  (i.J-1) 

-  A  f^  (1,J*1)  j  .  (16) 


4.  DISCUSSION 


The  supersonic  program  solves  accurately  the  partial  differential 
equations  derived  in  Section  2.  However,  we  should  note  that  the 
equations  (9a, b)  and  (10a, b)  are  approximate  representations  of  the  flow 
about  a  thin  wing.  Thus  the  numerical  results  should  be  viewed  as 
approximations  to  the  flow,  though  the  results  should  be  reasonably 
accurate,  provided  the  thickness  parameter  5,  is  small  and  the  Mach 
number  is  not  too  close  to  unity. 


The  program  offers  a  resolution  of  up  to  291  x  291  grid  points  on 
the  x-y  plane.  From  experience  with  the  program  5u  x  50  grid  points  gives 
excellent  agreement  (l.e.  four  significant  figures)  with  the  assymptotlc 
solutions  for  an  Infinite  rectangular  wing,  at  various  Mach  numbers  and  a 
sonic  triangular  wing. 


cfa; 


It  should  be  stated  that  the  program  uses  very  little  central 
processing  time.  When  the  functions  are  given  by  Table  1  and  the  data  by 
Table  2  (i.e.  the  sonic  triangular  wing)  the  running  time  is  7.2  c.p.u. 
seconds.  If  the  grid  step  length  is  halved,  that  is  put  DX  equal  to  0.01, 
then  the  running  time  is  30.6  c.p.u  seconds.  Also,  note  that  binding 
takes  about  six  c.p.u.  seconds. 


Finally  note  that  the  program  is  capable  of  handling  quite 
elaborate  wing  planforms.  For  instance:  forward  swept  wings  and  wings 
with  raonotonical ly  curved  or  convex  leading  and  trailing  edges;  however, 
the  program  will  in  most  oases  be  unable  to  determine  the  proper  grid 
system  for  wings  with  more  than  one  extreme  leading  or  trailing  point. 
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TABLE  1 :  Exatnple  of  function  subroutines 
FXL(lf),  FXT(Y).  HKX.Y)  and  H2(X,Y) 


FUNCT.F 


REAL  FUNCTION  FXL(Y) 

FXL=Y 

RETURN 

C  Alternate  entry  -for  -fxt. 

ENTRY  FXT(Y) 

FXT=1  .0 

RETURN 

END 

REAL  FUNCTION  HKX.Y) 
PARAMETER  ( PI=3 . 1 41 5926) 

DATA  THETA/0.5/ 

Hl=-1 .0»X*TAN(THETA*P1/180 ,0) 
RETURN 

C  Entry  -for  h2. 

ENTRY  H2<X,Y) 

H2=-l .0*TAN<THETA*PI/180 .0) 

RETURN 

END 


1 


TABLE  2:  Example  of  input  data  for  supersonic  program. 
The  first  two  parameters  have  Fortran  format  type  12, 
while,  the  rest  have  format  type  F12.8 


TESTl .D 


1 

ISW 

3 

I  OUT 

0.02 

OX 

1  .0 

B2 

1 .4142196 

AMS 

1 .4142196 

AME 

0 . 1 

DAM 

0 . 1 

AOMS 

0 . 1 

AOME 

0 . 1 

DAOM 

S' 

«- 

t 


APPENDIX  A 


SUPERSONIC. F 


PROGRAM  SUPERSONIC 

COMMON  /SUPER/  AL ,DX . B2,6S ,AM , IM , JC 
EXTERNAL  FXL , FXT . HI ,H2 

INTEGER  IS<29I ) , IE<291 ) , JSC 291 .291 > , JEC291 ,291 > 

DIMENSION  XLCSel ) .XT<58I > .XX(291 > 

COMPLEX  G(581 .146) ,F<291 ,291) ,FGC291 > ,ADUM 
CHARACTER«4  CHI C 4) ,CH2< 4) 

DATA  cm/' Phi  .'6"  ','Phi  ','Cp  "/ 

DATA  CH2/'Phi  , ’Q"  '  ,  ' Ph i ^ Cp-  '/ 

PARAMETER  <N=29I ,N2=58I ,N3=146) 

F0RMAT(2< 12/) ,8<F12.8/)) 

FO^AT  C '  1  '  ,5X , '  1 '  ,//  ,5X ,  'Unsteady  super  son  i  c  . '  ,//) 

FORMATC ' 1 ' ,5X  '2'  //, 5X Steady  supersonic . ' ,//) 

FORMATC ' 1 ' ,5X, '3' ,//,5X, 'Steady  supersonic  (symmetric).',//) 
FORMATC 1 IX, 'Grid  step  size  =  ' , El 0 . 5 . 1  OX , ' Semi  span  length  =  ', 
NF6. 4 ,// , 2X , ' Ini t iai  Mach  no.  =  ' ,F7.5,3X, 'Final  Mach  no.  =  ', 
♦*F7.5,3X,'Mach  no.  step  =  '  ,F7.5/) 

F0RMAT(2X, '  Ini  tial  -freq.  =  '  ,F7 . 5 , 2X , '  Fi  nal  -freq.  =  ',F7.5,2X, 
N'Freq.  step  =  ',F7.5/) 

FORMATC ' 1 ,///,6X, 13, 4X, II ,//,'  Mach  No.  =' , F7 . 5 , : , 5X , 
M'Frequency  =  ' ,F7.5) 

FORMATC//,'  y=  ' , F9 . 4 , 1  OX , ' Poi n ts  =' , I  4) 

FORMATC'  X  ' ,3X,10(Fa.5,4X)) 

FORMATC 1X,A4, IX, lOCEl 1 .4,1X)) 

Main  program  -for  supersonic  -flow  over  a  thin  airtoil. 

ISUI  =  1  -for  unsteady  supersonic  -flow 

ISW  =  2  -for  steady  supersonic  Flow  C  unsymmetr  i  c  ) 

ISU  =  3  -for  steady  supersonic  -flow  (symmetric) 

Output  para.  I  OUT  determines  type  o-f  output; 

I  OUT  =  1  -for  pressure  coef-ficient  only, 
lOUT  =  2  tor  potential  only, 

I  OUT  =  3  -for  both  potential  and  pressure. 

DX  =  step  size,  82  =  semi  span  length 

AMS  =  initial  Mach  no.,  AME  =  •final  Mach  no. 


STOP  '  Invalid  Mach  no.  range.' 
'  Invalid  switch  parameter  in 


DAM  =  del ta  Mach  no. 

I-t  ISW  =  1  then  AOMS  =  initial  ■freq.,  AOME  ==  ■final  ■freq. 
and  DAOM  =  delta  ■freq.. 

Leading  /trailing  edge  data  is  in  FXL/FXT . 

HI  represents  either  ho  or  d/dx.hs^f. 

H2  represents  either  d/dx.ho  or  d/dx.hs-. 

READ  C  5 , 4 1 , ERR=70 , END=7 1 ) I SW , I  OUT ,DX,B2,AMS,AME, DAM , AOMS , AOME , DAOM 

GO  TO  71 

CONTINUE 

STOP  '  Error  in  input  data.' 

CONTINUE 

IF  CAMS.LE.l .O.OR.AME.LE.l .0)  STOP  '  Invalid  Mach  no.  range.' 
IF  CISW.LT.l .0R.ISW.GT.3>  STOP  '  Invalid  switch  parameter  in 
input  data.' 

IF  CISU.EQ.l)  PRINT  42 

IF  CISW.EQ.2)  PRINT  43 

IF  CISU.EQ.3)  PRINT  44 

PRINT  45,DX,B2,AMS.AME,DAM 

IF  (ISW.EQ.l)  PRINT  46, AOiS, AOME, DAOM 

OUM=AME~AMS 

IF  CDUM.NE.0.0)  DAM-SIGNCDAM,DUM) 

AM-AMS 

CONTINUE 

AL=SQRT<AM»AM-I .0) 

CALL  GRIDXYCN,N2.XL,XT.IS,IE,JS,JE> 

IF  (ISW.EQ.l)  GO  TO  73 

AK-O.O 

AOM-0.0 

PRINT  47, JC, I  OUT, AM 
GO  TO  75 
CONTINUE 
DUM2-A0ME-A0MS 

IF  <DUM2.NE.0.0)  0A0M=-SIGN(DA0M,DUM2> 

AOM-iAOMS 

CONTINUE 

AK-Ajm#AM/AL/AL 

PRINT  47,JC,I0UT,AM,A0M 

CONTINUE  ’  ’ 


CALL  80LUE(I8U,N,N2,N3,JS,JE,AK.XL.XT,F,0> 
CALL  0UTPUTCN,t^2,XL,XT,XX,F,IS,IE,AK> 


OOOOOOO 


APPKNDIX  A  (coot.) 


SUPERSONIC. F 


IC=4-1ABS(5-2*ISU> 

IF  ( lOUT.EQ.l )  GO  TO  77 
00  76  J=1  ,JC 

IF  <IE<J>.LT.IS<J>>  GO  TO  76 
Y=FL0AT<2*J-1 )*DX/AL/2.0 
ITT=IE< J)-IS( J)+l 
NUM=INT<FLOAT<  ITT-t-8>/9 . 0  > 

IF  (NUM.EQ.O)  GO  TO  76 
PRINT  48, Y, ITT 
00  76  1*1 ,NUM 

JB=<l-n*9+lS<J> 
jT=MIN< IE< J) , 1*9+1 S<J)-1  > 

PRINT  49,<XX<K) ,K=JB, JT> 

PRINT  50  CH1(IC5,(REAL<F<K,  J)  >  ,t<=JB,  JT> 

IF  <ISW.NE.3)  PRINT  50 , CH2( I C) , ( A1MAG( F( K , J ) )  ,  K= JB  ,  JT  ) 
DNTINUE 
CON'.  INUE 

IF  <I0UT.EQ.2)  GO  TO  80 
IC=lC+t 

IF  (lOUT.NE.l)  PRINT  47 
AL..  ’=CMPLX<0 .0  ,-2.0»A0M*DX> 

00  J=1 ,JC 

IF  < IE( J) .LE. ISC J)+l .  GO  TO  79 

DO  78  I*IS< J>+1 .IE(J)-1 

FG<  I)*<F<  I-l  ,  J)-F(  I  +  l  ,  J>+AOUM*F‘’  I  ,  J>  )/DX 
CONTINUE 

Y=FLOAT< 2* J- 1 ) «DX/AL/2 . 0 
ITT=IE< J)-1S< J>-1 
NUM*INT<FLOAT< ITT+8)/9.0) 

PRINT  48, Y, ITT 
DO  79  1*1, NUM 

JB=<I-15*9+IS<J)+1 
JT=MIN<1E< J)-1,I»9+IS< J)) 

PRINT  49,<XX<KJ ,K=JB,JT> 

PRINT  50,CH1<IC>,<REAL<FG(K>),K=JB, JT) 

IF  <ISW.NE.3)  PRINT  50 , CH2< I C> , (AIlV^GC FG( K> > , K= JB , J 
CONTINUE 


GO  TO  74 


Mach  no.  step  size  is  zero.' 
iiND.DUM.NE.O.O)  GO  TO  72 


CONTINUE 

IF  (ISW.NE.l)  GO  TO  81 
AOM=AOM+DAOM 

IF  <DUM2*<AOM-AOME) .LE.O.O.AND.DUM2»DAOM.NE.O .0)  GOTO 
CONTINUE 

IF  COAM.EQ.O.O)  STOP  '  Mach  no.  step  size  is  zero.' 
AM=AM+DAM 

IF  <DUM*<AM-AME) .LE.O.O.AND.DUM.NE.O.O)  GO  TO  72 
STOP  '  Normal  completion.' 

END 

SUBROUTINE  6RI0XY<N,N2,XL,XT , I S, I E , JS , JE) 

CCMMON/SUPER/  AL,DX,B2,BS,AM.1M, JC 
INTEGER  IS<N) ,IE<N),JS(N,N>,JE<N,N) 

DIMENSION  XL<N2),XT<N2) 

Subroutine  to  initialize  grid  system. 

FIND  EXTREME  FORWARD  AND  TRAILING  POINTS 

Leading  and  trailing  edge  data  is  in  Function  subroutines 
FXL  and  FXT  <in  Eulerian  co-ordinates). 

DO  900  J-1.N2 
XL< J)-0.0 
XT<J)=0.0 
CONTINUE 
D2-DX/AL/2.0 
JJ8-1 
JJE-1 

e8*FXL<02) 

6E-FXT<02) 

XL(2)*GS 

>a<2)-GE 

J-2 

J-J+2 

Y-FL0AT<J-1)«D2 

IF  <Y.OT.B2>  00  TO  902 

XS-FXL<Y) 

XE*FXT<Y) 

XL< J>«X8 
XT<J)»XE 

IF  (XS.LT.OS)  THEN 


ooo  ono 
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GS=XS 
JJS=J/2 
END  IF 

IF  (XE.GT.GE)  THEN 
GE=XE 
JJE=J/2 
END  IF 
GO  TO  901 

902  CONTINUE 
JC=-l+J/2 

DO  903  J=2,2*JC.2 
XL< J)=XL< J)-GS 
XT< J)=XT< J)-GS 

903  CONTINUE 

DO  904  J*1,2*JC-1,2 
Y=FL0AT<  J-n*D2 
XL< J)=FXL<Y)-GS 
XT< J)=FXT<Y)-6S 

904  CONTINUE 
Y=FL0AT<2»JC)*D2 

IF  <Y.LE.B2)  THEN 

XL<2*JC-H  )=FXL<Y)-6S 
XT(2*JC+1 >=FXT<Y)-GS 
END  IF 

Find  leading  edge/characteristic  grid  points. 
IS< JJS)=1 

IF  (JJS.NE.l)  THEN 
DO  90<5 

ISC J)=IS< J+1 J+I 
DUM=XL  <2*J)/DX+2.0 

905  CONTINUE 

IF  <DUM.LE.FLOAT<IS( J>))  THEN 
ISC J)=IS<J)-1 
GO  TO  905 


f’06 


END  IF 
CONTINUE 
END  IF 

CJJS.NE.JC) 


IF 


THEN 


DO  908  J=JJS+1,JC 
ISC J)=ISCJ-l5+l 


DUM-XLC2»J)/DX+2.0 

907  CONTINUE 

IF  CDUM.LE.FLOATCISC J>)>  THEN 
ISC J)=ISC J)-l 
60  TO  907 
END  IF 

908  CONTINUE 
END  IF 

Find  trailing  edge/characteristic  grid  points 

lEC JJE)-1+INTCCGE-6S)/DX) 

IF  CJJE.NE.l)  THEN 

DO  910  J»JJE-1,1,-1 
IECJ)»IECJ+1)-I 
DUM»XTC2»J)/DX 

909  CONTINUE 

IF  CDUM.GE.FLOATCIECJ)))  THEN 
IECJ)-IE<J)+1 
60  TO  909 
END  IF 

910  CONTINUE 
END  IF 

IF  <JJE.NE.JC>  THEN 
DO  912  J-JJE4^1,JC 
IE<J)-IE<J-1>-1 
DUM-XT<2*J)/DX 

91 1  CONTINUE 

- - - 


IF 

<  DUM . OE . FLOATC I EC  J>  >  > 

IE<J)-IE<J)*1 

00  TO  911 

END 

IF 

912 

CONTINUE 

919 

END  IF 
J-JC+1 
J«J'l 

IF  (IE<J> 

.LT.ISCJ))  00  TO  919 

oooooo  OOOO  OOO 
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J.LT.l)  GO  TO  91.4 
OR.I .GT.IE< J>>  GO  TO  914 
60  TO  915 


JC=J 

FIND  TIP  DIAPRAGHM  REGION 

I J=INT<  < IE( JC)-IS< JC)>/2.0) 

JM=JC+IJ 

IF  < IJ.LE.O)  60  TO  918 
DO  913  J=JC+1,JM 
IS< J>=IS< J-1 )+l 
IE< J)  =  IE< J-1  >-l 

913  CONTINUE 

Find  j-startiriQ  and  stopping  grid  points  at  all 
ver  t i cal  1 evel s . 

918  IM=1E(JJE) 

DO  916  1=1 , IM 
JSUI=1 
J=0 

914  J=J+JSW 
IF  <  J , GT . JM . OR 
IF  <I.LT.IS<J) 

IF  <JSW.EQ.-1> 

JS< 1 , I )=2»J 
JSUI=-1 

J=JM+ 1 
GO  TO  914 

915  JE<l,n=2*J 

916  CONTINUE 

DO  917  K=2, IM 
K1=K-1 
IK=M0D<K,2) 

DO  917  1=1 ,IM-K+1 

JS<K, I )=JS<K1 , 1-IK+1)+2»1K-1 
JE<K, I )=MIN< JE<K1 ,1) ,JE<K1 ,I+1)>+1 

917  CONTINUE 
RETURN 
END 

SUBROUTINE  SOLVE< I SW ,N,N2 ,N3 , JS, JE , AK ,XL ,XT , F , 6> 

COMMON/SUPER/  AL . DX , B2,GS,AM, IM, JC 
INTEGER  JS<N,N) , JE<N,N5 
COMPLEX  F<N,N) ,G<N2,N3) ,DUM,DUM2 
DIMENSION  XL<N2),XT<N2) 

Subroutine  to  initialise  wing  boundari'  conditions,  grid  points 
and  integrate  the  unsteady/steady  supersonic  problem. 

Initialize  boundary  conditions. 

AOM-AK*AL*AL/AM 

AKM=AK»AM 

DD-DX/AL 

COEF=4.0/<8.0+AKeAK*DX»DX) 

DUM2«CMPLX<  0 .0,0.0) 

DO  498  I-1,IM 

XX-FLOAT < 2* I -3 ) *DX/2 . 0 
X-XX+GS 

DUM=CMPLX<COS<AKM*X) ,SIN<AKM*X> ) »DD 
JST=MAX< I . JS< 1 , I ) -1 ) 

DO  498  J-JST,JE<1 ,I)+1 ,2 
Y«FL0AT< J-1 ) eDD/2 . 0 
JJ-<J+l)/2 

IF  <Y.6T.B2.0R.XX.LT.XL<J>.0R.XX.QT.XT<J>>  THEN 
F<I ,JJ)»DUM2 

ELSE  IF  (ISW.EQ.l)  THEN 

F<I , JJ>-0UMeCMPLX<H2<X,Y> ,Hl<X,Y>*AaM> 

ELSE  IF  <ISU.EQ.2>  THEN 

F<I ,JJ)-0D*CMPLX(HI<X,Y>,-1 .0«H2<X,Y>> 

F<I , JJ)»ODeCMPLX(Hl<X,Y> ,0.0) 

END  IF 

498  CONTINUE 

DO  499  I-l ,IM 

DO  499  J-JS< 1,1) , JE< 1 , I ) ,2 
JJ* J/2 


499 


C0NrKi4'’'"‘'’“ 


A 
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C 

C 


550 


C 

C 

C 


Initialize  integration  points  to  zero. 

DO  550  1=1 ,N2 
DO  550  J=1 ,N3 

G(  I .  j:)=cmplx<o  .0 ,0 .0) 

CONTINUE 

Integrate  the  unsteady  and  steady  supersonic  problem. 
■1 


IST0P=2*IM- 
1 1=0 

551  11=11+1 
KE=IM-IAEtS(  I  I-IM) 
lK=i I I+l i/2 

DO  552  J=JS( 1 ,IK) ,JE<1 .IX) ,2 
J1=J+1 
J2=J-1 
JJ=J/2 

DUM=2.0«COEF*<6( J1 , 1 ) +G( J2 , 1 ) ) -G< J ,  1  >-COEF«F( IK,JJ) 

IF  < ISU.EQ.2.AND. ( J.GT.2*JC.0R.FL0AT( ( 1 1  + 1 >/2) . LT . 

1  XL<J)/DX))  THEN 

DUM2=0 ,5«<AIMAG<DUM>+REAL(DUM> ) 

DUM=CMPLX<DUM2,DUH2> 

END  IF 

F< IK, JJ)=DUM 
G< J, 1 )=DUM 

552  CONTINUE 

DO  553  K=3,KE,2 
IK=< I I-K+2)/2 
K1=<K+1 )/2 
K2=K1-1 

DO  553  «J=JS<K,IK),JE<K,1K),2 
J1=J+1 

G< J.Kl )=C0EF*<6< J1 ,K1)+G( J2,K1 >+G( J1 , K2> +G( J2 , K2>  > -6( J , K1 > 

553  CONTINUE 
11=11+1 

IF  < 1 1 .GT. ISTOP)  RETURN 
KE=1M-1ABS< 1 1-IM) 

DO  554  K=2,KE,2 
Kl=K/2 
K2=K1+1 
1K=< 1 l-K+2)/2 

DO  554  J*JS<K,1K) ,JE<K,IK> ,2 
Jl=J+l 
J2=J-1 

IF  < J • EQ . 1 )  J2=2 

G< J.Kl )=C0EF«<G< J1 ,K1)  +  G< J2,K1 >  +  G(Jl , K2> ♦Gt J2 ,K2>  > -G( J , K1  ) 

554  CONTINUE 
GO  TO  551 
END 

SUBROUTINE  0UTPUT<N,N2,XL .XT ,XX,F.IS,IE,AK) 

COMMON  /SUPER/  AL,DX,B2,GS,AM, IM, JC 
REAL  XL<N2) ,XT<N25 .XX<n5 
INTEGER  IS<N>,IE<NI 
COMPLEX  F<N,N5 
DO  930  1=1. IM 

XX<  I  )=FLd)AT<  I-l  )*DX+GS 

930  CONTINUE 
AKM— 1  .0*AK*AM 
DO  931  J=1 ,JC 

ARG-XL<2ij)/DX 

I1=2+INT<AR0> 

IF  <M0D<AR0.1 .0) .EQ.O.O)  11=11-1 

I2-INT<XT<2*J)/DX)+1 

IS<J>=I1 

IE(J>-I2 

DO  931  I-Il .12 

ARO=AKM«}OI(I)  _ _ 

F< I . J>-F< I , J>*CMPLX<COS<ARO> ,8IN» ARG) > 

931  CONTINUE 
RETURN 
END 
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SUBROUTINE  SUPDAT<M ,N , DX ,B2 , 1 SW ,  I  OUT , I C , AM , AOM , JC , JP , I  PI , 1 P2 , 
#Y1 ,Y2.X1 ,X2,F,G,*) 

DIMENSION  Y1 <M,N> ,Y2<M,N) ,X1 (M,N,N) ,X2<M ,N ,N> ,AM< M) ,AOM<M) 
INTEGER  IPl <M,N) , 1 P2(M ,N) , JC(M) ,JP<M) 

COMPLEX  F<M,N,N) ,G<M,N,N) 

F0RMAT<4X, I  1 ) 

F0RMAT<4</ ) ,28X,E10.5,29X,F6.4,9(/> , 7X , I  3 , 3X , I  1  , 1 2X , F7 . 5 , 

#17X,F7.5) 

FORMAT<4</) ,28X,E10 .5,29X,F6.4,7</) , 7X , I  3 , 3X , I  1  ,// , 1 2X , F7 . 5) 
FORMAT  <// ,  5)<; ,  F9 . 4 . 1  8X  !  1 4) 

F0RMAT<9X,9<F3.5,4X, :)) 

F0RMAT<6X,9<E1 1 .4,1X,:)> 

FORMAT < ) 

FORMAT<//> 

F0RMAT<7X, 13,//, 12X,F7.5,17X,F7.5) 

F0RMAT<7X, 13,//! 12X.F7.5> 

OPEN<  1  ,FILE=^SUPER.D'  , STATUS='' OLD'  ) 

IC=1 

READ<  1  .40  . ERR=29)  I SUI 

IF  (ISU.EQ.l)  READ< 1 ,41 ,ERR=29)DX,B2, JC< 1 ) , IOUT,AM< 1 ) ,AOM( 1 ) 
IF  <ISW.NE.l)  READ<1,42,ERR=29>DX,B2,JC<1> ,IOUT,AM(l) 
CONTINUE 

IF  (lOUT.EQ.n  GO  TO  24 
DO  23  J=1.JC<IC) 

READ<1 .43,ERR=29)Y2<IC,J> ,ITT 
NUM=INT<FL0AT<ITT+8)/9.0) 

IP2< IC, J)=ITT 
DO  23  K=1 ,NUM 
IS=K*9-8 
IE=MIN< ITT,K*9) 

REA0<1 ,44,ERR=29)<X2<IC,1,J) .I=1S,IE) 

READd ,45!eRR=29)<REAL(F<IC,I . J>) ,I=IS.IE> 

IF  <ISUI.NE.3)  READ<1,45,ERR=2^><AIMAG<F<IC,I ,J)),I=IS,IE) 
CONTINUE 

IF  <I0UT.EQ.2)  GO  TO  26 
READd  .46) 

CONTINUE 

DO  25  J=l,JCdC) 

READ<  1  ,43,END=28,ERR=29>Y1  dC,  J)  ,  ITT 
IF  dTT.llE.O)  GO  TO  27 
NUM=INT<FLOAT< ITT+8)/9.0) 

IPl  dC,  J)  =  ITT 
JPdC)=J 
DO  25  K=1 ,NUM 
IS=K*9-8 
IE=MIN< ITT.K*9) 

READCl  ,44,ERR=29)<XldC,I,J)  .I*IS,  IE) 

READ< 1 ,45!eRR=29) <REAL(G< IC, I . J) > , I=IS, IE) 

IF  dSW.NE.3)  READ<1  ,45,ERR=2^)<AlMAG<GdC,I  ,  J))  ,I  =  IS,  IE) 
CONTINUE 
CONTINUE 
READ<1,47) 

CONTINUE 

IC=IC+1 

IF  dC.QT.M)  GO  TO  30 

IF  <ISW.EQ.l)  READd  ,48,END=30,ERR=29)JC<  IC)  ,AMdC)  ,AOM(  IC) 
IF  (ISW.NE.l)  READd  ,49, END»30|ERR=29) JCdC)  JaMCIC) 

GO  TO  22 
CONTINUE 
IC=IC-1 
CONTINUE 
CL0SE< 1 > 

RETURN 
CONTINUE 
CLOSE < 1 > 

RETURN  1 
END 
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Description  of  Computer  Prograoi 


A  computer  program  based  on  the  finite  difference  scheme  (1^)  and 
(16)  has  been  Implemented  in  Fortran  on  the  ARL  ELXSI  computer.  The 
program  consists  of  a  main  program  called  SUPERSONIC  and  three 
subroutines,  GRIDXY,  SOLVE  and  OUTPUT.  For  completeness  a  listing  is 
supplied  in  Appendix  A.  The  position  _of  the  leading  and  trailing  edge, 
and,  either,  h^  and  hg^^  or  and  hg  must  be  supplied  by  the  user  in 
function  subroutines  called  FXL(Y),  FXTiY),  H1(X,Y)  and  H2(X,Y).  For  an 
example  see  Table  1  .  Note  that  X  and  Y  in  these  functions  are  sealed 
Eulerian  variables  and  not  Prandtl  -  Glauert  variables. 


The  program  has  three  options;  unsteady  flow,  steady  flow  and 

steady  symmetric  flow  (i.e.  hs  =  -  hs  x).  .  For  unsteady  flow 

Hl(x,y)  and  H2(x,y)  represent  hQ(x,y)  and  hQ^{x,y)  respectively.  The 
output  for  the  unsteady  case  is  the  real  and  imaginary  parts  of 

(XiY.o)  (see  (1(a))  and  Q(x,y,o),  where  the  unsteady  pressure 
coefficient  is  defined  as 


(P-  PJ 


Q  (x,y,o) 


iut 


(17) 


For  the  steady  case  H1  and  H2  represent  h  *x  and  h  x.  The  output 
for  this  case  is  s  s 


(x.y.o*),  *  (x,y,o  ),  (x,y,o*)  and  (x,y,o~), 

9  w  P  P 


where  In  the  steady  case. 


(P-  P.) 


Cp  (x.y.o*)  . 


(18) 
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For  the  steady  symmetric  case,  the  flow  below  the  x-y  plane  equals  the 
flow  above.^  The  output  for  this  case  is  Just  b,  (x,y,o  )  and 

C  (x.y.o  )  .  ^ 


After  FXL,  FXT,  HI  and  H2  have  been  programmed,  compile  tne 
functions  with  optimiser  switch  off,  that  is 

FORTRAN  MYFUNCTIONS  -OPT. 

Then  bind  the  programs  in  the  following  way: 


BIND  /USER/ST.  GEAR/SUPERSONIC , MYFUNCTIONS 

BFILE  =  MY  TITLE. 


The  executable  file  will  then  be  called  MYTITLE  and  it  will  exist  in  your 
area.  To  run  SUPERSONIC,  extra  data  must  be  supplied  through  an  input 
file.  Table  2  shows  a  typical  example  of  the  data  required  by  SUPERSONIC. 


The  first  parameter  in  Table  2  (ISW)  has  format  type  12  and  takes 
the  values  1,  2  or  3.  That  is  ISW  equals  one  for  unsteady  flow,  two  for 
steady  flow  and  three  for  steady  symmetric  flow.  The  second  parameter 
lOUT  also  has  format  type  12  and  it  also  takes  the  values  1,2  or  3.  It 
represents  the  type  of  output.  For  instance  if  lOUT  equals  one  then  only 
the  pressure  coefficient  is  output.  If  lOUT  equals  two  then  only  the 
potential  is  printed.  If  lOUT  is  three  then  both  the  potential  and 
pressure  coefficient  are  printed.  The  third  parameter  DX  has  format  type 


F12-8,  it  determines  the  distance  between  grid  points  on  the 
characteristic  grid.  Note  that  DX  is  twice  i  (3ee_(12a)).  When  FXL  and 
FXT  are  given  as  in  Figure  1,  the  Mach  number  is  /2  and  DX  -  0.02  as  in 
Table  2,  then  the  grid  system  on  the  x-y  plane  has  50  x  50  points.  Note 
that  SUPERSONIC  can  take  a  grid  size  on  the  x-y  plane  of  up  to  291  x  291 
points.  Also  note  that  if  an  "ACCESS  VIOLATION"  error  occurs  while 
running  SUPERSONIC,  then  the  grid  size  la  greater  than  291  x  291  points. 
To  remedy  this  situation  simply  increase  DX.  The  fourth  parameter  B2  also 
has  format  type  F12.8,  it  represents  the  sealed  Eulerian  semi-span 
length.  The  fifth,  sixth  and  seventh  parameters  determine  the  initial, 
final  and  delta  change  in  Mach  number.  That  is  AMS  is  the  initial  Mach 
number,  AME  la  the  final  Mach  number  and  DAM  Is  the  delta  change  in  Mach 
number.  Thus  with  a  suitable  choice  of  AMS,  AME  ai.d  DAM,  SUPERSONIC  will 
produce  data  for  several  different  Mach  numbers,  during  one  run.  Note 
that  AMS,  AME  and  DAM  have  format  type  F12.6.  Also  note  that  if  AMS  and 
AME  are  less  than  or  equal  to  one,  then  SUPERSONIC  will  not  execute.  The 
last  three  parameters  are  only  used  If  ISW  equals  1.  AOMS  represents  the 
Initial  reduced  frequency,  AONE  the  final  reduced  frequency  and  DAOM  the 
delta  change  in  reduced  frequency.  AOMS,  AOME  and  DAOM  have  format  type 
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F12.8.  Note  that  for-  each  Mach  number  AOMS,  AOME  and  DAOM  di.-to’ir.ine  a 
range  of  reduced  frequencies.  Also  note  that  if  ISW  is  not  equal  to  one 
then  these  parameters  need  not  be  supplied. 


Once  the  input  data  have  been  programmed  into  a  file  called,  for 
instance,  MYDATA,  then  SUPERSONIC  can  be  executed  by  entering: 


MYDATA  >  MYTITLE  >  MYOUTPUT, 


where  MYOUTPUT  represents  the  file  to  which  the  output  is  printed.  ,Note 
that  if  there  is  an  error  in  the  input  data  SUPERSONIC  will  stop  execution 
and  print  one  of  three  statements.  For  instance:  "invalid  switch 
parameter  in  input  data",  signifies  that  ISW  is  not  in  the  range  one  to 
three,  "invalid  Mach  no.  range",  signifies  that  AMS  or  AME  is  less  than  or 
equal  to  one,  or  "error  in  input  data",  signifies  an  error  condition  was 
encountered  during  the  input  operation.  If  SUPERSONIC  executes  normally 
then  upon  termination  SUPERSONIC  prints  "normal  completion"  to  your 
terminal  or  log  file. 


Table  3  shows  the  first  part  of  the  output  produced  when  the 
functions  of  Table  1  and  the  data  of  Table  2  are  used.  The  first 
character  on  the  first  line  is  a  Fortran  carriage  control  character  the 
second  character  is  ISW,  In  this  case  1.  DX  and  B2  are  printed  on  the 
sixth  line,  AMS,  AME  and  DAM  are  printed  on  the  eighth  line  and  if  ISW  -  1 
then  AOMS,  AOME  and  DAOM  are  printed  on  the  tenth  line.  On  the  twelfth 
line  (or  tenth  if  ISW  F  1)  the  first  character  is  a  Fortran  carriage 
control  character.  The  next  character  is  the  number  of  grid  points  in  the 
y  direction,  while  the  last  character  on  this  line  is  lOUT,  in  this  case 
3.  SUPERSONIC  then  prints  the  Mach  number  and  the  reduced  frequency  for 
the  first  run,  followed  by  the  value  of  y  at  the  grid  points  nearest  the 
wing  root  and  the  number  of  points  in  the  x-  direction  at  that  value  of  y, 
in  this  case  50.  SUPERSONIC  tjjen  prints  the  value  of  x,  and  the  real  and 
imaginary  parts  of  (or  ji  and  )  at  each  grid  point  for  that  y 
station.  After  all  sff  points  have  &en  printed  the  program  then  prints 
the  values  of  x  and  ^  at  the  next  value  of  y  until  the  wing  tip  is 
reached.  In  this  case  as  lOUT  equals  three  the  pressure  coefficient  at 
each  grid  point  (where  it  Is  calculated)  is  printed  In  the  same  way  as 
described  above.  Note  that  as  the  program  uses  central  differences  to 
calculate  the  x  derivative  of  ^  (or  the  pressure  coefficient  is 
not  calculated  on  the  leading  and 'trailing  edge  grid  points. 


To  facilitate  the  use  of  this  program  Appendix  B  contains  a  listing 
of  a  subroutine  that  can  be  used  to  read  the  output  of  SUPERSONIC.  When 
this  subroutine  is  used  the  number  of  configurations  read  is  transferred 
to  the  main  program  by  the  variable  IC.  The  Mach  number  and  reduced 
frequencies  are  stored  in  the  arrays  AM{K),  AOM(K)  where  1  S  K  S  IC. 


APPENDIX  C  (CONT.) 


The  number'  of  grid  points  in  the  spanwise  direction  for  the 
pressure/potential  is  stored  in  JP(K)/JC(K)  (1  S  K  S  IC).  The 
number  of  grid  points  in  the  streamwise  direction  at  the  Jth  spanwise 


point  for  the  pressure/potential  is  stored  in  IPI (K,J)/IP2(K,J) .  The 
value  of  y  at  the  Jth  spanwise  position  for  the  pressure/potential  is 
stored  in  If  1  (K , J)/if2(K,K) .  The  value  of  x  at  the  Jth  spanwise  and  Ith 
streamwise  position  for  the  pressure/potential  is  stored  in  Xi(K,I,J)  / 
X2(K,I,J)  and  finally  the  real  and  imaginary  parts  of  the  pressure 
coefficient  and  potential  at  the  Jth  spanwise  and  Ith  streamwise  grid 
point  are  stored  in  G(K,I,J)  and  F(K,I,J),  where  I  S  K  S  IC.  . 
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